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Abstract. This paper is devoted to the analysis of a numerical scheme for the coagulation and 
fragmentation equation with diffusion in space. A finite volume scheme is developed, based on 
bJT)' a conservative formulation of the space nonhomogeneous coagulation-fragmentation model, it is 

shown that the scheme preserves positivity, total volume and global steady states. Finally, several 
. numerical simulations are performed to investigate the long time behavior of the solution. 
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Coagulation and fragmentation processes occur in the dynamics of cluster growth and describe 
the mechanisms by which clusters can coalesce to form larger clusters or fragment into smaller 
i— —i | ones. Models of this type have found important applications in areas ranging over polymer kinetics 

|29j . aerosols |16| . cluster formation in astrophysics [26] to the animal grouping in biology |27j . 
This paper is devoted to the numerical simulation of the dynamics of the density function 
(NJ ' / = f(t,x,y) > of particles (polymers, clusters) with position x G Q C 1^ (N > 1), volume 

y G R + := (0,oo), and time t > 0. The distribution function / is subject to coagulation and 
fragmentation phenomena with respect to the volume variable y and diffusion in space x G ft and 
is governed by the following equation 

(1) % ~ d(y)A x f =Q(f), 

O 

where Q(f) is the coagulation-fragmentation operator acting on the volume variable y G M + and 
these particles diffuse in the environment f2, which is assumed to be a smooth bounded domain 
with normalized volume |f2| = 1. Moreover, the model ([I]) is supplemented with an initial datum 

(2) f(0,x,y) :=r(x,y), 

and in the following analysis we will consider homogeneous Neumann boundary conditions on d£l 

(3) V x f{t,x,y) ■ u(x) = 0, (t,x,y) G R + x <90 x R+ 

with v the outward unit normal to S7. We also assume the diffusion coefficient d(y) to be non 
degenerate in the sense that there exists d* ,d* G M + such that 

(4) d* > d(y) > 4 > 0, y G R + . 

Let us first describe precisely the coagulation and fragmentation process. In the simplest coagulation- 
fragmentation models the clusters are usually assumed to be fully identified by their size (or vol- 
ume, or number of particles). The coagulation- fragmentation models we consider in this paper 
describe the time evolution of the cluster size distribution as the system of clusters undergoes 



binary coagulation and binary fragmentation events. More precisely, denoting by C y the clusters 
of size y with y G M + , the basic reactions taken into account herein are 

(5) C y + Cy/ ^'^ C y+y i, (binary coagulation) 
and 

(6) C y ^ UiV) C y -yt + Cyt, (binary fragmentation) , 

where a and b denote the coagulation and fragmentation rates respectively, and are assumed to 
depend only on the size of the clusters involved in these reactions. 
Thus, the coagulation-fragmentation operator is defined by 

(7) Q(f) = Q c (/)-Qf(/), 

where the coagulation operator is 

I rv r°° 

(8) Qc(f)(x,y) = - / a{y',y - y') f{x,y') f(x,y -y')dy' - / a{y,y') f{x,y) f{x,y') dy' 

z Jo Jo 

whereas fragmentation mechanism by which a single particle splits into two pieces is given by 

i r-y poo 

(9) Qf(f)(x,y) = - / b{y' : y-y')dy' f(x,y) - / b(y, y') f(x, y + y') dy' . 

1 Jo Jo 

The coagulation coefficient, a = a(y,y'), characterizes the rate at which the coalescence of two 
particles with respective volumes y and y' produces a particle of volume y + y', whereas the 
fragmentation coefficient, b = b(y,y'), represents the rate at which the fragmentation of one 
particle with volume y + y' produces two particles of volume y and y' . Both coefficients a and b 
are nonnegative functions and 

a(y, y') = a(y', y), b{y, y') = b(y' , y), 



(10) 



a, b G Lf oc (K + x R+) . 



For symmetric kernels, we observe that during the microscopic coagulation and fragmentation 
processes, as depicted in equations ©-(H]), the number of particles varies with time while the 
total volume of particles is conserved. 

In terms of /, the total number of particles and the total volume of particles at time t > are 
respectively given by 



M (t,x):= [ f(t,x,y)dy, M x (t,x) := f yf(t,x,y)dy. 

JR+ JR+ 



Besides existence and uniqueness results, few is known on the qualitative behavior of solutions 
of coagulation-fragmentation equations except when coagulation and fragmentation coefficients 
are linked by the so called "detailed balanced condition" : there exists a nonnegative function 
M eL^fix R+, (1 + y)dxdy) such that 

(11) a(y,y')M(y)M(y') = b( y ,y')M(y + y'), (x, y, y') G O x M+ x M+. 

This condition implies that M is a stationary solution to the coagulation-fragmentation equation 

d(y)A x f + Q(f) = 0. 
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We further assume that, for each R € M + , the equilibrium M satisfies the positivity condition 



(12) inf M(y) > 0. 

ye[0,R] 

An additional and interesting feature of the detailed balanced condition is that there exists an 
entropy H given by 



H(f\M) 



dy dx, 



'QxM+ 

which satisfies the following H theorem 

^ = —DU) - [ { ^P^(t,x,y)dydx < 0, 
where the entropy dissipation of the coagulation-fragmentation operator is defined by 
D(f) =11 («//'- bf") (ln(a//') -ln(6/'0) dydy'dx, 

with the shorthand notation / := f(t,x,y), f := f(t,x,y') and /" := f(t,x,y + y'). Since D(f) 
only vanishes when / is an equilibrium, it naturally means that 

f(t) — ► M, when t -> oo. 

The main purpose of this work is to present a numerical scheme to solve (pQ) built upon a finite 
volume discretization with respect to the space variable x £ £1 and volume variable y G M + . The 
analysis of the so-obtained scheme would allow to prove the convergence of the discretized particle 
density towards a solution to the continuous problem on a fixed time interval [0, T] (T > 0), see for 
instance [1], but here our aim is different. It will consist in the study of the long time behavior of 
the numerical solution on a fixed mesh in space and volume variables. Indeed, often in applications 
we are interested in steady states or in the long time behavior of the solution, it is then important 
to design a numerical scheme, which has good stability properties uniformly in time and remains 
consistent with respect to the exact solution for long time. Thus, the scheme proposed in this 
paper is designed such that it preserves qualitative properties of the exact solution as steady 
states of the coagulation-fragmentation operators (|8])-([9]). Moreover, an estimate of the entropy 
dissipation due to an appropriate finite volume approximation with respect to space and volume 
variables is given and the scheme is shown to give a consistent approximation of the continuous 
problem in the long time asymptotic limit t — > oo. 

Before describing more precisely our results, let us recall that the coagulation and fragmentation 
equations (pQ) with ([8|)-(|9]) have been the object of several studies recently. 

On the one hand, among the various approaches for the approximation of coagulation and 
fragmentation models, we may distinguish between deterministic and Monte Carlo methods. We 
refer for instance to [U [T2l [23] for deterministic methods, [21 [17] for stochastic methods, and the 
references therein. Concerning the convergence analysis of numerical methods for coagulation 
and fragmentation models, we refer to [22] for a rigorous study of quasi Monte-Carlo methods. 
For deterministic approximations, the situation is different since the relationship between discrete 
and continuous models has been considered by some authors, see the survey paper [8] and pQ. 
A rigorous setting for the formal analysis performed in PQ under general assumptions on the 
coagulation and fragmentation coefficients has been given in [21]. Then, similar techniques are 
used to prove convergence of discrete schemes to the exact solution in [1] . 
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However, few is known concerning the stability and the analysis of the numerical solution in 
the long time asymptotic limit. For the continuous model, in reference [20], the authors are able 
to apply the techniques of weak compactness, as in the paper by R.J. DiPerna and P.-L. Lions [7] 
about the Boltzmann equation, to prove the weak stability of weak solutions. In the case when 
an .H-theorem holds, they can obtain some partial information about the large-time behavior of 
solutions. More recently, J. A. Carrillo et al. prove exponential decay of the solution towards 
equilibrium for the inhomogeneous Aizenman-Bak model [5] using entropy dissipation methods 
[6]. Our analysis for discrete models will be inspired by these works. 

We now briefly outline the contents of the paper. In the next section, we introduce the numerical 
approximation of ([1]) and state the stability result which we prove in Section [3] In the final section 
(Section [5|), some numerical simulations are performed with the numerical scheme presented in 
Section [2] and the long time behavior of the solution is investigated. 



2. Numerical scheme and main results 



In order to compute an approximation of this model using a finite volume method in space 
variable x E O, and volume variable y G R + , we observe that the coagulation- fragmentation 
operator can be written in a conservative form. Indeed, writing equation (pQ) in a "conservative" 
form, as proposed in [261 [28] , enables to describe precisely the time evolution of the total volume. 
Also, this formulation is particularly well adapted to a finite volume discretization which, in 
turn, is expected to give a precise account of volume conservation. Precisely, the coagulation and 
fragmentation terms can be written in divergence form: 



yQc(f)(x,y) 

yQt{f){x,y) 
where the operator C(f) is given by 



dy 
dy 



x,y), 



(13) C(f)(x,y) : 
and J^if) is 

(14) Hf )(*,!/) 



u a(u, v) f(x, u) f(x, v) dvdu , (x, y) G 0, x 



y-u 



u b(u, v) f(x, u + v) dvdu , (x, y) G fl x 



Finally, the coagulation-fragmentation equation reads 

df „ , , , dC(f) 



(15) 



V 



d(y)yA x f 



+ 



8t dy 
f(0,x,y) = f in (x,y), (x,y)en,xl 



dJFU) 
dy 



and we assume that the initial datum / in is a nonnegative function which satisfies: 
(16) / in G L x (n x R+) HL 1 ^ x R + ,ydxdy). 

Here and below, the notation L l {Q, x M + , ydxdy) stands for the space of the Lebesgue measurable 
real-valued functions on x M + which are integrable with respect to the measure ydxdy. 
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When designing the volume discretization of the coagulation and fragmentation terms, it is 
necessary to truncate the infinite integrals in formulae (|13j) - (|14p . But this means restricting the 
domain of action of kernels a and b to a bounded set of volumes y, that is, preventing coagulation 
to occur among particles with volume exceeding a fixed value. The discretization we propose is 
based on the conservative truncation method for the coagulation and fragmentation terms. Given 
a positive real R, let (x, y) E f2 X (0, R) 



C R (f)(x,y) 



R-u 



u a(u,v) f(x,u) f(x,v) dvdu . 



(17) 



x—u 
R 



u a(u, w — u) f(x, u) f(x, w — u) dw du . 



In that case, C R (f)(x,0) = C R (f)(x, R) = foreach x G SI so that the total volume of the 
solution is now nonincreasing with respect to time. 

As regards the fragmentation term, the truncation is also a conservative truncation on the 
fragmentation term. Using the same idea, we introduce for (x, y) £ $7 x (0, R) 

"R-u 

u b(u, v) f(x, u + v) dv du. 

iy-u 

ry t-R 

(18) = / / ub(u,w — u) f(x,w) dv du. 

Jo Jy 

Then, the conservative coagulation-fragmentation operator satisfies exactly the conservation of 
total volume, so that the following equation is indeed a truncated conservative coagulation and 
fragmentation equation: 

dC R (f R ) 



(19) 



since 



V 



dfR 
dt 

f(0,x,y) 



d{y) y A x /r 
= f in (x,y), 



d_ 

dt 



dy + 
(x,y) enx (o,R), 



dF R {f R ) 
dy 



(x,y) 



V fR(t,x,y)dxdy = 0. 



Under the detailed balance condition (fTTI) . model (TL9j) has also a steady state Mr, only depending 
on y and such that 

M R (y) = M(y), ye[0,R}. 

and 

■L I I r{x,y)dxdy= [ M R {y)dy. 

\ lL \ Jo Jn Jo 
Convergence for large values of R has been thoroughly studied in the recent past. We briefly 
mention some results for the coagulation equation (that is, with 6 = 0). These results adapt easily 
to the coagulation-fragmentation equation but under different assumptions on the kernels. Let us 
mention that when a(y, y 1 ) /{y y') — > as y + y' — > +oo, convergence as R — > +oo of the solutions 
to (fT9|) toward a solution of (fT3|) - (fT5|) can be proven by using the approach developed in [21]. 
Thus, since the convergence of solutions to (fT9|) towards solutions of (fT3]) - (fl5|) is well established 
in rather general situations, this paper will only focus on the convergence of a sequence built on 
a numerical scheme towards a solution to the equation (|19|) when the truncature R is fixed. In 
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the remainder of the paper, for the sake of clarity, we drop the subscript R and write / instead 
of fn for a solution of equation (|19p . Parameter R being fixed, this should raise no confusion. 

Now, we turn to the discretization of equation (|19p. Having reduced the computation to a 
bounded interval, the second step is to introduce the space and volume discretizations. To this 
end, we set Ay G (0, 1) and N y a large integer, and denote by (l/i-i/2)ie{o,...,N } a mesh of (0, R), 
where = i&V, and A* = [yi-i/ 2 , yi+1/2) for i > 0. 

Concerning the space variable, we choose an admissible mesh of given by a family T of control 
volumes (open and convex polygons in 2-D), a family £ of edges and a family of points (xk)k&t 
which satisfy Definition 5.1 in [llj. It implies that the straight line between two neighboring 
centers of cells (xk,xl) is orthogonal to the edge a = K\L. In the set of edges £, we distinguish 
the interior edges a G £% n t and the boundary edges a G £ e xt- For a control volume K G T, we 
denote by £k the set of its edges, £% n t,K the set of its interior edges, £ e xt,K the set of edges of K 
included in V = d£l. 

In the sequel, we denote by d the distance in M. N , m the measure in M. N . We assume that the 
family of mesh considered satisfies the following regularity constraint : there exists £ > such 
that 

(20) d(x K ,a) > £d(x K ,x L ), for K £ T, for a e£ int ,K, a = K\L. 
The size of the mesh is defined by 

(21) 5 = max(diam(K)) . 

For all a E £ , we define the transmissibility coefficient: 

m(cr) 



d(xR,XL) 

m(er) 



for a G £ int , a = K\L, 
for a £ £ ex t,K- 



Next, X(T) will be the set of functions from f2 to R which are constant over each control volume 
K G T. 

We define the approximation f h (0) of the initial datum / m as usual by 

Ny 

(22) f h (o) = EE/few,, 

i=0 KeT 

with 



A,; 



where 1e denotes the characteristic function of the subset E of O and converges strongly to f° 
in L 1 ^ x (0, R)) as /i = (Ay, 5) goes to 0. 

Let us now introduce the numerical scheme itself. For each integer i G {0, • • • , N y } and each 
K G T, we define the approximation of f(t,x,y) for i G M + and (x,y) G if x Aj as fx,i(t)- 
The sequence {fK,i)K,i is defined by the following discretization of the coagulation-fragmentation 



equation : for K G T, i £ {0, . . . , -/V^}, we solve the ordinary differential system 

Vk, 
dt 



(23) 



We have set 



m (^) ~n~ ~ d (Vi-l/2) ^2 T ° D K,afK,i = m(K) Q Kji , 



fK.m := ff. 



( 24 j = r 1 t 

Vi-i/2^y y%-i/2&v 

where the flux Cj(,i+i/2 in (f2i|) represents an approximation of the coagulation operator fT7|) . 
J r K,i+i/2 is the approximation of the fragmentation part (|18p and both are defined by 

i N y -l 

(25) C K,i+i/2 = ^2 ^2 A V 2 Vj-i/2 a j,i-j Ikj fK,l-j, 

j=0 l=i+l 

where a^j := a{yi,yj) with yi = {i + 1/2) Ay and 

i Ny-l 

(26) ^K,i+\/2 = ^2 ^2 A y 2 yj-i/2 b j,i-j fK,i, 

j=0 l=i+l 

where bij := b(yi,yj) whereas the fluxes at the boundary are 

(27) Cft-,-1/2 = F K _i/ 2 = C KtNy+1 / 2 = ^K,N y +i/2 = 0) K e T. 

Concerning the approximation of the diffusion in space, we set vk,<t the unit normal to a outward 
from K and define an approximation of V x / • v^ a on a by 



(28) D K><T f K>i 



fh,i — fx,i, if cr = K\L G £i n t,K, 
0, if a G £ ex t,K, 



for all K G T . This discretization obviously relies on a simple finite volume approach for the 
space and volume variables (see, e.g. [4j). 

The following function f h defined on M + x ft x [0, R] will be useful in the sequel. 

N v 

(29) f\t,x,y) = J2 £/tf,i(*)lxxA,- 

A'eT i=0 

We may now state our main result. 

Theorem 2.1. Assume that the coagulation and fragmentation kernels satisfy H]), H0\) and (112)) 
and / in satisfies H6\) . We consider a uniform volume mesh in y and require the mesh T in space 
to satisfy condition if20j) . (TJip. 

Then, there exists a unique solution f h to the finite volume scheme ( 123)) - (128 )) . which satisfies 



(i) nonnegativity 



f K At)>0, (K, i) G T x {0, . . . , N y }, te 
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(ii) conservation of total volume 



d_ 

dt 



m(K) Ay ?/j_ 1/2 fK,i(t) 

K.i 



0, V* 6 



(iii) entropy dissipation 
dH(f h \M) 



dt 



where 



and 



+ ^ E E d (yi-i/2)T*D K>a f K>i ln(^ = -D{f h ), 

K,i <t&£ k ^ 1 ' 

fx/ 



H(f h \M) = ^m(K)Ay 



A',, 



fK,i In 



1 +Mi 



a M-z fx,i Ik. 



h,i-l fx,i 



D (f h ) = ^EE A ^ m ( K ) (%i-lfK,lfK,i-l ~ In 
Z=0 ■ 

(iv) Moreover, the finite volume scheme is asymptotic preserving i.e. 

f h (t) — ► M h , ast — ► oo. 

where 

M h (x,y) = Mi := exp(-ayi), V(x,y) G if X J/i+1/2) 
and a G M + is such that 

]T ^m^i-l^ = J] ^m(K)y,_ 1/2 /^. 
AeT i=0 AeT i=0 



3. A PRIORI ESTIMATES 

As we mentioned in the introduction, our goal is not to prove that the sequence of functions 
(f h )h converges in some sense to a function / as h = (Ay, 5) goes to 0, but to study the long 
time behavior of the numerical solution and to prove its uniform stability with respect to time. 
First, we prove that the solution f h to the scheme ([23p -(|27 p enjoys properties similar to those of 
function / given by (|19|) which we gather in Proposition 13.21 below. 

On the one hand, we first re-write the discrete coagulation-fragmentation obtained form the 
conservative form (|19p in a new form which is consistent with the classical formulation dTJ-Q. 
On the other hand, we find a discrete version of the weak formulation from which we prove a 
priori estimates. 



Proposition 3.1. Assume the approximation of the coagulation-fragmentation model 1 11)-^W\I 
is given by the finite volume scheme li23\) -(23\ ) . Then, the discrete operator satisfies 



Q K ,i 



2/1-1/2 A v 



+ 



F A,i+l/2 — F. A,i-l/2 

Vi-1/2 &y 



with 

1 ' 

(30) Q K ,i = - E A V ( a 3,i-3 fK,j SK,%-3 ~ bj,i-j I 



3=0 

Ny-1 

— E Ay (dij-i fx,i fK,j-i - bi,j-ifK,j) 
j=i 

Proof: Starting from the definition of the fluxes (|25p and (|26|) , we set 

~ Ck,i+1/2 ~ CK,i-l/2 F K,i+l/2 ~ ? . K,i-X/2 

UK,i = T 1 T 

Vi-l/2 A y Vi-1/2^V 

and by construction 

i Ny-l 

Vi-i/2QK,i = ~E E A yyj-i/2 ( a j,i-j fi<j fi<,i-j - fx,i) 

3=0 l=i+l 

+ EE A yVj-l/2 (Hl-j fKj fKJL-j ~ b j,l-jfK,l), 
3=0 l=i 

which can be easily simplified into 

Ny-1 

Z/i-1/2 Q-K,i = - E A yyi~l/2 ( a i,l-i fK,i fK,l-i - t>i,l-ifK,l) 
l=i+l 
i-1 

+ E Ay yi~w* ( a hi-i fx,i fx,i-i - h,i-i fK,i) ■ 

1=0 

Then, adding and subtracting the term j/i_i/ 2 (<H,o fx,o fx,o - h,o fK,i), we get 

i 

!/t-l/2 Qlf.i = ( a l,i-l Ik,1 fK,i-l — k,i-lfK,i) 

1=0 

Ny-1 

— E Ay yi-l/2 (. a i,l-i fK,i fK,l-i - h,l-ifK,l)- 
l=i 

Hence, it only remains to treat the first term of the right hand side: we compute for a symmetric 
sequence g^j such that g^j = g 3 \i for each G {0, • • • ,N y } 2 

i i i 

}Jgi,i-i = /jgi,i-i - y~](i-l)gi,i-i 

1=0 1=0 1=0 

i i 

= E* 5i '*- j ~ E Zft -M- 

z=o z=o 
Thus, using the symmetry g^_i \ = gi j_; we finally get 



j i i 

1=0 1=0 
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Therefore, with y^in = i Ay and taking successively g^j = djj fK,ifK,j & n d gtj = bij, it yields 
the result 



1 1 

2=0 



l=i 

□ 



Next, we prove the following estimates 



Proposition 3.2. Assume the approximation of the coagulation-fragmentation operator (17)- 
U8\) is given by the finite volume scheme jl25\)- (2B\) . Then, for each K E T the discrete operator 
(Q-K,i)o<i<N y satisfies a discrete weak formulation : for any sequence (<Pi)o<i<N v an d f or eac -h 

K eT 



Nv-1 



(31) ^ AyQ K)i (pi = -- ^2 ^2 Ay 2 (ai,i-i fx,i fK,i-i - h,i-i fx,i) (<?z + ¥>i-; - ¥>i) 

8=0 i=0 Z=0 

Proof: We multiply ([30j) by Ay c/?j and sum over i G {0, ■ ■ ■ , iVy — 1} to get 

(32) Ayg^ m = Ay 2 (Jj - J 2 ) , 



i=0 



with 



i=0 j=0 

Ny — X Ny — 1 

h = / J ( a i,j-i ffC,i fK,j-i - bij-if'K,j) <fi- 

i=0 j=i 

We keep the first term I\ as it is and only treat the term li- First, we commute indexes (i, j) and 
then perform a change of index on I2 which yields 

Ny-i j 

h = S ^{ a i,j-i fK,i fK,j-i — bi,j-ifK,j) ^Pi 

j=0 i=0 

Ny-l i 

(33) = /if,i-z /r:,z — h-1,1 fic,i) fi-i- 

i=0 «=0 

Using the symmetry Ojj = o^j and 6jj = 6^ for each (i, j) £ {0, • • • ,N y — 1}, we also have that 

Ny-l i 

(34) I 2 = ^2 /, i a l,i-l fK,l fK,i-l ~ h,i-lfK,i) <Pi- 

8=0 2=0 
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Finally, gathering (f3"2"|) - (f34"|) . we get the result 



• 1 Ny-l i 

^2 AyQ K!i ipi = -- ^ ^2 A V 2 ( a l,i-l fx,l fx,i-i - h,i-l fK,i) (<Pi + Vi-i ~ <Pi 

i=0 i=0 1=0 

□ 

Now we are ready to prove Theorem 12.11 

4. Proof of Theorem 12.11 

We assume that the kinetic and diffusion coefficients fulfil (|4]), (jlOl) and (112j) . respectively, and 
that a and b are positive a.e. in R + x M + . We are also given an initial datum _f in satisfying (I16p 
and denote by f h the solution to (i23"j) constructed in (f2"4"j) - (j26j) . 

On the one hand, existence and uniqueness of a solution to 

(35) - d(yi_i/ 2 ) ^ T a D K ,afK,i = m(K)Q K ,i, 

o-e£K 

directly follows by applying the classical Cauchy-Lipschitz theorem for a finite set of ordinary 
differential equations. Existence of a global solution will be a consequence of the following uniform 
estimates with respect to time. 

On the other hand, since the discrete distribution function is solution to (|35|) . the nonnega- 
tivity of f h (t) follows from the monotonicity of the discrete operator ^<re£ T(J ^K,afK,i and the 
structure of the discrete coagulation-fragmentation operator (|30p written as the sum of a gain 
operator and a local loss term. 

Next, we prove the total volume conservation 

, Ny-l 

(36) - £ 52m(K)Ay yi __ 1/2 f K)i (t) = 0. 

at i=0 KeT 

Indeed; we multiply ([33]) by Ay J/i_i/2 and sum over (K, i) G T x {0, . . . , N y — 1} 

^^2 m ( K ) ^VVi-l/2fK,i{t) = d(yi-l/ 2 ) ^ T v D K,ofK,iDK,o{yi-l/2) 

K,i K,i itG£x 

+ ^2 m ( K ) Ay QkM-i/2- 

K,i 

Thus using a discrete integration by part, we prove that the first term of the right hand side is 
zero, whereas applying Proposition 13.21 with cp^ = yi-1/2, it also gives that 

^m(K)Ay Q K ,m-i/2 = 0, 

K,i 

which concludes the proof of (ii). 

Now, let us prove the stabilization towards an equilibrium in the long time when the kinetic 
coefficients a and b satisfy the detailed balance condition (jlip . As mentioned in the introduction 
the detailed balance condition (jlip ensures an analogue of the Boltzmann //-theorem for the 
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coagulation-fragmentation equations which we derive now for the discrete solution to (|23|) . We 
set 

H(f\M) = Y,m(K)Ay (f K>i (hx(^) - l) + M, 



KA 



and take in the discrete weak formulation (|3ip . the test function <p such as ipi = \a{fK,i/Mi) with 
Mi = M(y i _i/ 2 ) and noticing that 



(f h ) = Y,m(K)Ay 



KA 



dt 



In 



Mi 



and recalling pip with (pi = hi(fjc,i/Mi), we obtain 

dH 
~dt 



(f h ) + Ay E ^-1/2)^^/^111 



-Ay J2 m ( K ) &K,i In 



/gj 



which yields 
(37) 



"eft" 



TS ■ r-C JK\LA 

K,i a£tn 1 ' 



with fK\L,i = Uka - fL,i)/Qn(fK,i/ fL,i)) > 0. Then, we treat the right hand side using the 
detailed balance condition (llip and get 



D (f k ) = ^EE A ^ Ki-2 fK,l fKA-l - \i-ljK, 



KA 1=0 



hi 



dlA-l flQ, fKA-l 
bli-l fK,i 



which proves (Hi). 

Now, we study the long time asymptotic behavior of the numerical solution f h . To this aim 
we establish the following estimates 

Lemma 4.1. The numerical solution f h satisfies the following uniform estimates with respect to 
time t £ R + : there exists kq > such that 

52m(K)Ay(l + yi _ 1/2 )f Kti (t) + H(f h (t)\M) < k q . 

KA 

Moreover, for all t G R + , we have 



(38) 
and 



Ay 



E E d(Vi-i/2) 1/2 H<r)\DfK,U*) 

i=0 crG£ int 



ds < Kq. 



(39) Ay 2 ^EE ( a iA-i fK,i fKA 



-I — bli-l / 



KA 



KA 1=0 



111 



O-lA-l flQ, fKA-l 
h,i-l fKA 



ds < Kq, 
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Proof: Using the total volume conservation and the entropy dissipation, we have already proven 
that 

^MK)Ayy^ 1/2 f K 4t) + H(f h (t)\M) < C . 

K,i 

It remains to show that there exists a constant Co > 0, only depending on the initial datum f m 
and M, such that 



(40) 



^m{K) Ay f Kti (t) < Co- 



K,i 



K,i 



Indeed, on the one hand we notice that 

Y,m(K) Ay f K>i (t) < e^m^Ai/l^^^M,}^ 

K,i 

+ \ E m ( K ) A V 1 {/K, ! W>e 2 Mi} Sk,x In 

K,i \ i / 

< e 2 m ( K ) A V M i 

K,i 

+ ^m(K)A, fc ln( M 

is i \ 



On the other hand, observing that r ln(r) > r | ln(r)| — 2/e for r > 0, we have 



fK,i In 



Mi 



In 



Mi 



2 Mi 



hence it gives 



Ayf K , 

K,i 



fK,i(t) 

Mi 



f K At) 

Mi 



< Y,m(K)Ayf Kii ln 

K,i 

e ^— ' 

< H (f h \M) + m(K) Ay (f Kti + * M^j . 

Combining the two inequalities yields 

^mWAyfc(t) < (e 2 + ^)^(K)AyM i + \H(f h (t)\M), 
K,i ^ ' K,i 

^ ( e2 + -) E m W A ^ M * + \H(f h (0)\M) =:C . 



K,i 
13 



Next, let us prove ([39]) . We start with the entropy estimates ([37|) . which we integrate with 
respect to time t G M + , it gives 

(DK, a fK,iY 



H{f h {t)\M)+ / A^E E d fa-i/2K 

< ff(/*(0)|M). 



+ / D(f h (s))ds 

JK\L,i JO 



Then, since H is a convex function, we know that H(f h (t)\M) is bounded from below and we 
get the result (p9|) . 

Finally, we show there exists Cq >, only depending on / m and M, such that 

2 



(41) 



i=0 o-e^int 

(T~K \ L 



ds < C , Vi G 



To this aim, we start with the Cauchy-Schwarz inequality 

Ny 

E E M°)d(Vi-i/2) 1/2 \DfKte\, 



i=Q o&£i n t 
a = K\L 



< 



f Ny \Df K - \ 2 \ 1/2 ( Nv 

E E ^%i-i/2) /'^ E E m ( a 

i=0 <re£i„i J-ftT|L,! I I j=0 CTS£int 

\ a = K\L / \ cr=K\L 



1/2 



r)d(x K ,X L ) f K \ L>i 



Applying ([20]) and ([36]), ([3ZD, it follows that 
/ N v 



i=0 o-6£ int 
a=K\L 



2 

< - 



E E Aym^Xy^) 1 / 2 !^, 
=o 

/ 

V 



\ 2 

/ 



E E w( yi _ 1/2 )B^ z 

a=K\L / N 



E E *vMK)fK,i ] , 

if eT 



< 



2 Cp 



E E A y T ° d (yi-i/2. 



\Df K , 



V 



= CT££ int 

<T = if |i 



Then, for each i G M + we integrate the latter inequality with respect to time over the interval 
[0, t], and use the entropy dissipation (137D which guarantees that the right hand side is uniformly 
bounded with respect to time. □ 



Now, we prove that the numerical solution converges to a steady state. Let (t n ) n >o be a 
sequence of positive real numbers such that t n — ► +oo and set f n (t) = f h (t n + 1) for n > 1 and 
t G M + . Owing to the construction of f n , it is easily seen that f n is a weak solution to (|23p - ([2l 
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on [0, +00) with initial datum j (t n ). We fix T G M + and infer from above that 

sup Vm(iO Ay (l + R) + H(f n (t)\M) < k . 
te[o,T] K i 

Moreover, by construction of f n and applying Lemma 14.11 it implies that 
(42) 



A y / E E rf&^i/a^m^iD^coicft 

< Ay / t+T V V ^-v^rmVHAfo >il<r (t)\dt 



and 
(43) 



«=0 aeSi n t 



0. 



n — > 00 



A ^ 2 / EE fh fh-i - b ^ fh) 

J ° K,i 1=0 
r-t n +T i 

< Ay 2 / EE ^ ai ' i - 1 f K < 1 ~ hl ^~ l ? K 



111 



b l,i-lfx,i )_ 

%t-l fK,l fK,i-l 
h,i-l fK,i 



dt 



dt 



K,i 1=0 

— ► 0, as n — > 00. 

Thanks to theses estimates, we deduce that there are a subsequence of (/ n )neN (not relabeled) 
and a weak solution / to (|2"3"l) - (|27l) such that 

(44) f n - /, mC([0,T]), asn -» +00. 
Moreover, passing to the limit in (|30p . it easily follows that 

(45) QkAH — > Qjf,i(/), m^^T), asn^+oo. 
and passing to the limit in the discrete diffusion operator, we get 

(46) d{yi_ 1/2 ) ^2 T ^ D fK,i,a — * d(yi_ 1/2 ) ^ T a Dfx,i,a in L 1 (0, T) , asn^+cx). 

On the one hand, we use (f4*3j) and (|44|) to conclude that 

£>(/) = 0, o.e.in (0,T). 

Consequently, / satisfies 

a i,j fK,i(t) fK,j(t) = bi,j fK,i+j{t), 

for almost every t G (0, T) and each (K, G T x {0, . . . , -^V y } 2 . In particular, it implies that 

Gjf,i(/(t)) = 0, GTx{0,...,iVj 2 . 
On the other hand, (f42l) and (l46j) ensure that 

(47) = 0, (K, i,a)eTx{0,...,N y }x£ 

for almost every t G (0, T). Therefore, /(i) does not depend on time and satisfies the steady 
states equation (j23|) . Moreover, using the zero flux boundary conditions (|3|) and ([4"7]) . / does 
not depend on K G T. Recalling that the differential equation ()23[) conserves global volume, we 
conclude there is a G M + such that 

] K ,% ■= Mi, a.e. in (0,1% 
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where 

Ny-l Ny-1 

i=0 K£T i=0 KeT 

5. Numerical simulations 

This section is devoted to the numerical study of the convergence to equilibrium under the 
detailed balance condition and when this condition is not satisfied. We also investigate numerically 
the case with non homogeneous Dirichlet boundary conditions. 

5.1. Detailed balance kernels and convergence to equilibrium. We assume that the co- 
agulation and fragmentation coefficients fulfil the detailed balance condition: there exists a non- 
negative function M G Lj;(IR + x R + ), such that 

(48) a(y,y')M(y)M(y > ) = b(y,y')M(y + y% (y, y') G R + x R+. 

Since in that case there exists a Lyapunov functional H at the discrete level, we have proven that 
the numerical solution f h converges to a discrete equilibrium M h . Here, we choose kernels a and 
b as follows: 

a(y,y') = Hy,y') = i, 

which yields M(y) = exp(— yj^/M\). It is the Aizenman-Bak model for reacting polymers which 
diffuse in space with a non degenerate size-dependent coefficient d(y) = a > 0. In [5], the authors 
demonstrate that the entropy-entropy dissipation methods developed by Desvillettes and Villani 
for the Boltzmann equation [6] applies directly and gives the exponential convergence with explicit 
rates towards global equilibrium for constant diffusion coefficient in any spatial dimension or for 
the non degenerate diffusion in dimension one. Thuerefore, the global equilibrium is given by 

M(y) = expi-y/y/Mj, y G K + , 
where the number Mi is given by 

Ml := / / yf in (*,y)dxdy, 

and M satisfies 

d(y) A X .M + Q(M) = 0. 
On the other hand, the local equilibrium Mi oc is 

Mi oc (x,y) = exp(-y/ v / Mi(x)) 

where the function x G £1 — > M\(x) G R + is given by 

Mi(x) := / yf(t,x,y)dy, 
ii+ 

and Mi oc satisfies 

Q(M loc ) = 0. 

The relative entropy H{f\M) can be split in two different parts 

(49) H(f\M) = H{f\M loc ) + H(M loc \M), 

where the first term in the right-hand side represents the "distance" between / and the local 
equilibrium whereas the second term evaluates the distance between the local and global equilibria 

16 




and only depends on the macroscopic quantity M\. In the following we present some numerical 
results. As an initial datum, we take 

f in (x,y) = exp(-a(x)y), 

with a(x) = 1 + 0.1 cos(2 ir x\) cos(2 it x%) and x = {x\, X2) G (—1/2, 1/2) 2 . The diffusion is taken 
to be constant d(y) = 0.1, R = 20, N y = 64 and next 128 and finally At = 0.002. 

Let us first mention that in [14\ 115] a similar problem i.e.; trend to equilibrium of the solution 
to the nonhomogeneous Boltzmann equation, is investigated numerically. It is shown that the 
relative entropy with respect to the local equilibrium oscillates with time when the solution 
becomes close to the equilibrium. Here, we will show that the situation is completely different 
and much simpler. 

In Figure [H we report the evolution of the total number of particles Mq, the high order mo- 
ments in y of /" with respect to time and the Lyapunov functional H(f\M). As expected, the 
total volume M\{t) remains constant throughout time evolution and the moments stabilize to a 
fixed value. As regards the asymptotic profile, our numerical results are in fair agreement with 
the equilibrium M(x) = exp(— x). Moreover, we observe that the scheme is able to give the 
correct behavior of the numerical entropy H(f\M), which converges exponentially fast to zero [5]. 
However, we observe that trend to equilibrium for the coagulation-fragmentation with diffusion is 
much more simpler than trend to equilibrium for the Boltzmann equation [T3] since no oscillation 
occurs for the relative entropy H(f\Mi oc ). Indeed, after a short transient regime, the solution / 
converges to equilibrium as an exponential with respect to time. 



5.2. Convergence to equilibrium for not detailed balance kernels. After this first result, 
we now explore a different situation where nothing is known about equilibrium and entropy 
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(a) (b) 

Figure 2. Evolution of (a) the total number of particles Mq, the total volume Mi, 
M2 and M3 (b) the functional H{f\M), H(f\Mi oc ) and f h (t ~ 00, y) in log scale. 



dissipation. Indeed, we choose kernels a and b as follows: 

a(y,y') = (yy) 1/2 , b{y,y') = 1. 

The situation becomes much more complicated since we do not know neither the expression of 
the steady state nor the expression of entropy. As an initial datum, we take 

P n {x,y) = exp(-a(x)y), 

with a(x) = l + 0.5cos(47rxi) cos(47T3;2) and x = (xi,x%) G [— 1/2, 1/2] 2 . Moreover, the diffusion 
is taken to be 

d(y) = 0.1/(1 + y), 

which is degenerated for large y. We choose the truncation such that R = 20, N y = 64 and 
next 128 and finally At = 0.002. In Figure El we still report the evolution of the total number 
of particles Mq and other high order moments in y of f h , that is Mi, M2 and M3. As expected, 
the total volume M\{t) remains constant throughout time evolution, whereas high order moments 
vary and next stabilize to a fixed value, which means that the solution converges to a steady state. 
In the same figure (right hand side), we also report the distribution function with respect to y in 
log scale for large time, which indicates that the tail with respect to y of the steady state / is 
still exponentially decreasing and / is of course constant in x G 

5.3. Convergence to equilibrium for non homogeneous boundary conditions. In this 
last section, we study the coagulation-fragmentation operator with diffusion in space and mixed 
boundary conditions in x G More precisely, we aim to approximate the solution in = 
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(0,1/8) x (0,1) to 

' ^ - d(y)A x f = Q(/)„ 
fit = 0, x, y) = / ln (x, y), (x, y) G ft x R+, 

< 

/(i, xi = 0, x 2 , y) = exp (-y/a(x 2 )) , (x 2 , y) G <9ft x R+, 

V x / • u(x) =0, xi = 1/8, or X2 = 0, or x 2 = 1, 

with a(x 2 ) = (1 + cos(47rx 2 ))/2 and v the external unit normal to <9ft. 

We consider an initial condition, which is at equilibrium for the coagulation fragmentation 
operator 

/ in (x,y) = exp(-y/a(x)), (x,y) G ft x R+, 
with a(x) = (1 + cos(32 7rxi) cos(47rx 2 ))/2. 

We perform numerical simulations using the finite volume scheme with diffusion coefficient 
d(y) = 0.01/(1 + y) and a = b = 1. The following figures (Fig [3], H]and[5]) are illustrations of the 
profile of the number of particles Mo(i, x) and total volume Mi(t,x) for t G M + and x G ft and 
also the projections 

/•1/8 

P(t,x 2 ,y)= yf(t,x,y)dx u x 2 G (0, 1), y G K+. 

JO 

It allows us to observe the qualitative behavior of the system. The profiles are in that case smooth 
and also converges to an equilibrium. For such a computations we have considered N y = 64 and 
128 x 128 grid points for the space variable x G ft. In this situation, we cannot characterize 
explicitly the equilibrium in space since the diffusion coefficient is not constant and H is not 
valid here due to the Dirichlet boundary conditions for x\ = and then the total volume is not 
preserved at all. However, we still observe that the numerical solution converges to a discrete 
equilibrium (see Fig. [3]i5]). 

6. Conclusion 

In this paper, we make use of different principles (volume conservation, entropy dissipation, 
existence of steady states) which allow to build a stable and accurate numerical scheme for the 
nonlinear dynamics of coagulation, fragmentation and diffusion equations. Such an algorithm is 
able to recover the main properties of the exact solution and in some particular cases, it is proven 
that the method is asymptotically stable, in the sense that the numerical solution converges to an 
equilibrium which is consistent with the exact equilibrium of the system. Numerical simulations 
illustrate the efficiency of the algorithm even when we cannot prove convergence to equilibrium. 

Here, we have only considered the time evolution of a distribution function, but the method 
can be easily coupled with Euler or Navier-Stokes equations after some adaptations. Typical 
applications are transport problems (including linear and nonlinear diffusion or fluid dynamics), 
which are coupled with population balance dynamics represented by a distribution function / 
depending on space x G ft and "size" variable y G M + . 

On the other hand, for some coagulation and fragmentation kernels [12} [4"]. total volume is not 
conserved at all and then a non conservative formulation can be used to dicretize the coagulation 
and fragmentation operators using the same kind of finite volume scheme. Then, the finite volume 
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Figure 3. Time evolution of the density Mq{x) at time t = 0; 0.33; 0.66 and 4 

approach allows to use non uniform meshes for the volume variable y € R + , which is particularly 
well suited in this case |12j . 
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